{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\newcommand{norm}[1]{||#1||}$\n",
    "$\\newcommand{mbx}{\\mathbf x}$\n",
    "$\\newcommand{mbu}{\\mathbf u}$\n",
    "$\\newcommand{d}{\\text{div}}$\n",
    "$\\newcommand{mbf}{\\mathbf f}$\n",
    "$\\newcommand{t}{\\text{tr }}$\n",
    "$\\newcommand{mcA}{\\mathcal A}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 线弹性问题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模型"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\mu\\Delta\\mbu - (\\lambda+\\mu)\\nabla\\d\\mbu = \\mbf\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A(\\mbu) = \\frac{\\nabla\\mbu + \\nabla\\mbu^T}{2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\sigma(\\mbu) = 2\\mu A(\\mbu) + \\lambda(\\t A(\\mbu)) I\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\mcA(\\sigma) = \\frac{1}{2\\mu}\\left(\\sigma - \\frac{\\lambda}{n\\lambda + 2\\mu}(\\t\\sigma) I\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\d\\sigma(\\mbu) = \\mbf\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 存在问题"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 解法器效率太低\n",
    "* 矩阵组装基本是线性复杂度，但仍有改进空间\n",
    "* 需要更多测试例子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 测试模型数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display\n",
    "from sympy import *\n",
    "import numpy as np\n",
    "init_printing()\n",
    "\n",
    "def linear_elasticity_model(u, v):\n",
    "    n = len(v)\n",
    "    d, e = symbols('lam, mu')\n",
    "    du = u.jacobian(v)\n",
    "    Au = (du + du.transpose())/2\n",
    "    sigma = d*Au.trace()*eye(n)+2*e*Au\n",
    "    f = -diff(sigma[:, 0], v[0])\n",
    "    for i in range(1, n):\n",
    "        f -= diff(sigma[:, i], v[i])\n",
    "    return f, du, sigma, d, e\n",
    "\n",
    "def simplify_linear_elasticity_model(u, v):\n",
    "    n = len(v)\n",
    "    du = u.jacobian(v)\n",
    "    sigma = du + du.transpose()\n",
    "    f = -diff(sigma[:, 0], v[0])\n",
    "    for i in range(1, n):\n",
    "        f -= diff(sigma[:, i], v[i])\n",
    "    return f, du, sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例 1** 求解区域为 $\\Omega = [0, 1]^3$\n",
    "$$\n",
    "\\mbu = \\begin{pmatrix}\n",
    "2^4 \\\\ 2^5 \\\\ 2^6\n",
    "\\end{pmatrix}\n",
    "x(1-x)y(1-y)z(1-z)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[-a*x*y*z*(-y + 1)*(-z + 1) + a*y*z*(-x + 1)*(-y + 1)*(-z + 1), -a*x*y*z*(-x + 1)*(-z + 1) + a*x*z*(-x + 1)*(-y + 1)*(-z + 1), -a*x*y*z*(-x + 1)*(-y + 1) + a*x*y*(-x + 1)*(-y + 1)*(-z + 1)], [-b*x*y*z*(-y + 1)*(-z + 1) + b*y*z*(-x + 1)*(-y + 1)*(-z + 1), -b*x*y*z*(-x + 1)*(-z + 1) + b*x*z*(-x + 1)*(-y + 1)*(-z + 1), -b*x*y*z*(-x + 1)*(-y + 1) + b*x*y*(-x + 1)*(-y + 1)*(-z + 1)], [-c*x*y*z*(-y + 1)*(-z + 1) + c*y*z*(-x + 1)*(-y + 1)*(-z + 1), -c*x*y*z*(-x + 1)*(-z + 1) + c*x*z*(-x + 1)*(-y + 1)*(-z + 1), -c*x*y*z*(-x + 1)*(-y + 1) + c*x*y*(-x + 1)*(-y + 1)*(-z + 1)]])\n"
     ]
    }
   ],
   "source": [
    "x, y, z, a, b, c= symbols('x, y, z, a, b, c')\n",
    "u0 = x*(1-x)*y*(1-y)*z*(1-z)\n",
    "u = Matrix([a, b, c])*u0\n",
    "f = linear_elasticity_model(u, (x, y, z))\n",
    "print(f[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct linear system time: 1.659951514000113\n",
    "Solve time: 0.08165926999981821\n",
    "Construct linear system time: 13.403858605000096\n",
    "Solve time: 2.987838785999884\n",
    "Construct linear system time: 108.33645204000004\n",
    "Solve time: 192.936152407\n",
    "Ndof: [  750  4374 29478]\n",
    "error: [[  1.39739596e+00   2.29570086e-01   2.73310407e-02]\n",
    " [  1.46080383e+01   4.92832365e+00   1.25911061e+00]\n",
    " [  9.33096445e-02   8.14716287e-03   5.08085642e-04]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**例 2** 求解区域为 $\\Omega = [0, 1]^3$\n",
    "\n",
    "$$\n",
    "u = \\begin{pmatrix}\n",
    "e^{x-y}x(1-x)y(1-y) \\\\\n",
    "\\sin(\\pi x)\\sin(\\pi y)\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[-lam*(x*y*(-x + 1)*(-y + 1)*exp(x - y) - 2*x*y*(-y + 1)*exp(x - y) + 2*y*(-x + 1)*(-y + 1)*exp(x - y) - 2*y*(-y + 1)*exp(x - y) + pi**2*cos(pi*x)*cos(pi*y)) - 2*mu*(x*y*(-x + 1)*(-y + 1)*exp(x - y) - 2*x*y*(-y + 1)*exp(x - y) + 2*y*(-x + 1)*(-y + 1)*exp(x - y) - 2*y*(-y + 1)*exp(x - y)) - 2*mu*(x*y*(-x + 1)*(-y + 1)*exp(x - y)/2 + x*y*(-x + 1)*exp(x - y) - x*(-x + 1)*(-y + 1)*exp(x - y) - x*(-x + 1)*exp(x - y) + pi**2*cos(pi*x)*cos(pi*y)/2)], [-lam*(-x*y*(-x + 1)*(-y + 1)*exp(x - y) - x*y*(-x + 1)*exp(x - y) + x*y*(-y + 1)*exp(x - y) + x*y*exp(x - y) + x*(-x + 1)*(-y + 1)*exp(x - y) - x*(-y + 1)*exp(x - y) - y*(-x + 1)*(-y + 1)*exp(x - y) - y*(-x + 1)*exp(x - y) + (-x + 1)*(-y + 1)*exp(x - y) - pi**2*sin(pi*x)*sin(pi*y)) - 2*mu*(-x*y*(-x + 1)*(-y + 1)*exp(x - y)/2 - x*y*(-x + 1)*exp(x - y)/2 + x*y*(-y + 1)*exp(x - y)/2 + x*y*exp(x - y)/2 + x*(-x + 1)*(-y + 1)*exp(x - y)/2 - x*(-y + 1)*exp(x - y)/2 - y*(-x + 1)*(-y + 1)*exp(x - y)/2 - y*(-x + 1)*exp(x - y)/2 + (-x + 1)*(-y + 1)*exp(x - y)/2 - pi**2*sin(pi*x)*sin(pi*y)/2) + 2*pi**2*mu*sin(pi*x)*sin(pi*y)]]) \n",
      "\n",
      "Matrix([[x*y*(-x + 1)*(-y + 1)*exp(x - y) - x*y*(-y + 1)*exp(x - y) + y*(-x + 1)*(-y + 1)*exp(x - y), -x*y*(-x + 1)*(-y + 1)*exp(x - y) - x*y*(-x + 1)*exp(x - y) + x*(-x + 1)*(-y + 1)*exp(x - y)], [pi*sin(pi*y)*cos(pi*x), pi*sin(pi*x)*cos(pi*y)]])\n"
     ]
    }
   ],
   "source": [
    "x, y = symbols('x, y')\n",
    "u = Matrix([exp(x-y)*x*(1-x)*y*(1-y), sin(pi*x)*sin(pi*y)])\n",
    "f, du, sigma= linear_elasticity_model(u, (x,y))\n",
    "ss, cc, e, t0 = symbols('ss, cc, e, t0')\n",
    "sub = {sin(pi*x)*sin(pi*y):ss, cos(pi*x)*cos(pi*y):cc, exp(x-y):e, (-x+1)*(-y+1):t0}\n",
    "print(f, '\\n')\n",
    "print(du)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**例 3** 求解区域 $\\Omega = [0, 1]^2$\n",
    "\n",
    "$$\n",
    "u = \\frac{\\pi}{2}\\begin{pmatrix}\n",
    "\\sin^2(\\pi x)\\sin(2\\pi y)\\\\\n",
    "-\\sin^2(\\pi y)\\sin(2\\pi x)\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[-lam*(-pi**3*sin(pi*x)**2*sin(2*pi*y) - 2*pi**3*sin(pi*y)*cos(2*pi*x)*cos(pi*y) + pi**3*sin(2*pi*y)*cos(pi*x)**2) - 2*mu*(-pi**3*sin(pi*x)**2*sin(2*pi*y) - pi**3*sin(pi*y)*cos(2*pi*x)*cos(pi*y)) + 2*pi**3*mu*sin(pi*x)**2*sin(2*pi*y) - 2*pi**3*mu*sin(2*pi*y)*cos(pi*x)**2], [-lam*(2*pi**3*sin(pi*x)*cos(pi*x)*cos(2*pi*y) + pi**3*sin(2*pi*x)*sin(pi*y)**2 - pi**3*sin(2*pi*x)*cos(pi*y)**2) - 2*mu*(pi**3*sin(pi*x)*cos(pi*x)*cos(2*pi*y) + pi**3*sin(2*pi*x)*sin(pi*y)**2) - 2*pi**3*mu*sin(2*pi*x)*sin(pi*y)**2 + 2*pi**3*mu*sin(2*pi*x)*cos(pi*y)**2]]) \n",
      "\n",
      "Matrix([[pi**2*sin(pi*x)*sin(2*pi*y)*cos(pi*x), pi**2*sin(pi*x)**2*cos(2*pi*y)], [-pi**2*sin(pi*y)**2*cos(2*pi*x), -pi**2*sin(2*pi*x)*sin(pi*y)*cos(pi*y)]])\n"
     ]
    }
   ],
   "source": [
    "x, y = symbols('x, y')\n",
    "u = pi/2*Matrix([sin(pi*x)**2*sin(2*pi*y), -sin(pi*y)**2*sin(2*pi*x)])\n",
    "f, du, sigma, d, e= linear_elasticity_model(u, (x,y))\n",
    "print(f, '\\n')\n",
    "print(du)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 简化模型求解\n",
    "\n",
    "$$\n",
    "-\\Delta \\mbu - \\nabla \\d \\mbu= \\mbf\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\sigma (\\mbu) = \\nabla\\mbu + \\nabla\\mbu^T\n",
    "$$\n",
    "\n",
    "$$\n",
    "-\\d\\sigma (\\mbu) = \\mbf\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**例 2** 求解区域为 $\\Omega = [0, 1]^3$\n",
    "\n",
    "$$\n",
    "u = \\begin{pmatrix}\n",
    "e^{x-y}x(1-x)y(1-y) \\\\\n",
    "\\sin(\\pi x)\\sin(\\pi y)\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Matrix([[-3*x*y*(-x + 1)*(-y + 1)*exp(x - y) - 2*x*y*(-x + 1)*exp(x - y) + 4*x*y*(-y + 1)*exp(x - y) + 2*x*(-x + 1)*(-y + 1)*exp(x - y) + 2*x*(-x + 1)*exp(x - y) - 4*y*(-x + 1)*(-y + 1)*exp(x - y) + 4*y*(-y + 1)*exp(x - y) - pi**2*cos(pi*x)*cos(pi*y)], [x*y*(-x + 1)*(-y + 1)*exp(x - y) + x*y*(-x + 1)*exp(x - y) - x*y*(-y + 1)*exp(x - y) - x*y*exp(x - y) - x*(-x + 1)*(-y + 1)*exp(x - y) + x*(-y + 1)*exp(x - y) + y*(-x + 1)*(-y + 1)*exp(x - y) + y*(-x + 1)*exp(x - y) - (-x + 1)*(-y + 1)*exp(x - y) + 3*pi**2*sin(pi*x)*sin(pi*y)]]) \n",
      "\n",
      "Matrix([[-pi**2*cc - 3*e*t0*x*y + 2*e*t0*x - 4*e*t0*y - 2*e*x*y*(-x + 1) + 4*e*x*y*(-y + 1) + 2*e*x*(-x + 1) + 4*e*y*(-y + 1)], [e*t0*x*y - e*t0*x + e*t0*y - e*t0 + e*x*y*(-x + 1) - e*x*y*(-y + 1) - e*x*y + e*x*(-y + 1) + e*y*(-x + 1) + 3*pi**2*ss]]) \n",
      "\n",
      "Matrix([[x*y*(-x + 1)*(-y + 1)*exp(x - y) - x*y*(-y + 1)*exp(x - y) + y*(-x + 1)*(-y + 1)*exp(x - y), -x*y*(-x + 1)*(-y + 1)*exp(x - y) - x*y*(-x + 1)*exp(x - y) + x*(-x + 1)*(-y + 1)*exp(x - y)], [pi*sin(pi*y)*cos(pi*x), pi*sin(pi*x)*cos(pi*y)]])\n"
     ]
    }
   ],
   "source": [
    "x, y = symbols('x, y')\n",
    "u = Matrix([exp(x-y)*x*(1-x)*y*(1-y), sin(pi*x)*sin(pi*y)])\n",
    "f, du, sigma= simplify_linear_elasticity_model(u, (x,y))\n",
    "ss, cc, e, t0 = symbols('ss, cc, e, t0')\n",
    "sub = {sin(pi*x)*sin(pi*y):ss, cos(pi*x)*cos(pi*y):cc, exp(x-y):e, (-x+1)*(-y+1):t0}\n",
    "print(f, '\\n')\n",
    "print(f.subs(sub), '\\n')\n",
    "print(du)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 快速算法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "$$\n",
    "\\begin{pmatrix}\n",
    "A & B^T\\\\ B & -C\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "x_0 \\\\ x_1\n",
    "\\end{pmatrix}\n",
    "= \\begin{pmatrix}\n",
    "0 \\\\ b\n",
    "\\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{pmatrix}\n",
    "A & B^T\\\\ B & -C\n",
    "\\end{pmatrix}^{-1}\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "I & A^{-1}B^T \\\\ 0 & -I\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "A & 0 \\\\ B & S\n",
    "\\end{pmatrix}^{-1}\n",
    "$$\n",
    "\n",
    "其中 $S = BA^{-1}B^T+C$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{pmatrix}\n",
    "D & B^T\\\\ B & -C\n",
    "\\end{pmatrix}^{-1}\n",
    "= \n",
    "\\begin{pmatrix}\n",
    "I & D^{-1}B^T \\\\ 0 & -I\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "D & 0 \\\\ B & \\tilde S\n",
    "\\end{pmatrix}^{-1}\n",
    "$$\n",
    "\n",
    "其中 $\\tilde S = BD^{-1}B^T+C$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "核心需要实现的计算是：\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "D & 0 \\\\ B & \\tilde S\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "u_0 \\\\ u_1\n",
    "\\end{pmatrix}\n",
    "=\\begin{pmatrix}\n",
    "r_0 \\\\ r_1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "易知\n",
    "\n",
    "$$\n",
    "u_0 = D^{-1} r_0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "进而用代数多重网格可求解下面方程：\n",
    "\n",
    "$$\n",
    "\\tilde S u_1 = r_1 - BD^{-1}r_0:=r_2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 1:** GS 前磨光\n",
    "$$\n",
    "u_1^1 = u_1^0 + G( r_2 - \\tilde S u_1^0)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 2：** 粗空间 AMG 校正\n",
    "$$\n",
    "u_1^2 = u_1^1 + \\Pi M \\Pi^T(r_2 - \\tilde S u_1^1)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Step 3：** GS 后磨光\n",
    "$$\n",
    "u_1^3 = u_1^2 + G^T(r_2 - \\tilde S u_1^2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
